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USING  SIMULATION  TO  ESTIMATE  FIRST  PASSAGE  DISTRIBUTIONS 

by 

Sheldon  M.  Ross  and  Zvi  Schechner 

0.  INTRODUCTION 

Consider  a  discrete  time  Markov  process  {Xr  ,  n  ■  0,1,  ...}  such  that 
whenever  the  present  state  is  x  the  next  state  is  chosen  according  to  the 
distribution  Px  .  Let  Xq  ■  0  be  fixed  and  consider  for  a  given  set  of 
states  A  the  number  of  transitions  N  until  the  Markov  process  enters 
the  set  A  .  We  are  interested  in  estimating  the  distribution  and  the  mean 
of  N  by  use  of  simulation. 

Assuming  that  the  process  {Xr}  is  such  that  N  is  finite  with  Prob¬ 
ability  1,  this  model  can  be  simulated  by  the  standard  Monte  Carlo  method  of 

letting  Xq  -  0  and  then  generating  a  random  variable  from  the  distribu¬ 

tion  Pq  to  determine  X^  .  If  X^  ■  x  ,  we  then  generate  a  random  variable 
from  the  distribution  Px  and  set  it  equal  to  ^  ,  and  so  on.  We  stop  the 
run  when  we  obtain  a  state  in  A  .  The  process  is  then  repeated  for  a 
second  run,  and  so  on  until  a  total  of  r  (a  fixed  predetermined  number) 
simulation  runs  are  obtained.  If  we  let  N^  ,  i  -  1,  ...,  r  denote  the 
number  of  steps  until  A  is  reached  in  the  ith  run,  then  the  usual  way  to 
estimate  P(N  ■  k)  ,  k  >_  1  and  E[N]  is  as  follows: 

P(N  •  k)  is  estimated  by  P{1  :  N^  «  k}/r 

E[N]  is  estimated  by  £  N^/r  . 

i-1 

In  this  paper,  we  propose  new  estimators  for  these  quantities  which  are 
based  on  the  "observed  hazards"  also  known  as  the  "predictable  projection". 


That  is,  consider  a  given  simulation  run  terminating  at  N  and  define 


y 


where  is  the  (n  -  l)th  state  and  Px(A)  is  the  probability  that  the 

next  state  visited  from  x  is  in  A  .  We  use  the  process  {Xn>  to  es¬ 
timate  the  various  quantities  of  interest.  In  Section  1,  we  present  two  es¬ 
timators  for  E[N]  .  The  first  of  these  can  be  used  to  form  an  interval 
estimate  and  we  show  that  it  is  unbiased,  having  a  smaller  variance  than 
the  classical  Monte  Carlo  estimator.  The  second  estimator  presented  is  a 
point  estimate  of  E[N]  .  In  Section  2,  we  consider  the  problem  of  estimating 
the  distribution  of  N  .  In  Section  3,  we  consider  the  distribution  of  the 
final  state  X^  .  Section  4  deals  with  the  continuous  time  analog  and  as  a 
consequence  a  new  formula  for  estimating  convolution  is  obtained. 
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ESTIMATING  THE  MEAN  NUMBER  OF  TRANSITIONS 


1.1  Method  1 


For  a  given  simulation  run  if  we  let  N  denote  the  number  of  transi¬ 
tions  needed  to  enter  A  ,  then  N  can  be  expressed  as 


»  “  I  h 


where 


if  j  £  N 
if  j  >  N 


Hence, 


E[N]  -  l  E[I  ]  =  1  +  l  EflJ 

J-l  J  j-2  3 


and  the  usual  estimator  can  be  regarded  as  estimating  E [ Xj  J  by  looking  at 
the  values  of  1^  for  that  cycle  and  then  summing  over  j  to  estimate  E[N] 
We  propose  to  estimate  1^  by  conditioning  on  the  first  j-2  transi¬ 


tions  .  As 


I  X0,Xr  Xj_2]  “  >_  j  |  Xq, 


1  -  PY  (A)  if  N  >  j  -  2 
j-2 

0  if  N  £  j  -  2 


Hence,  our  estimate  for  E[N]  based  on  a  single  simulation  run  ending  when 
the  process  enters  A  is 


"l"  L"  .  ■  I  ' 


N+l  r  *1  N+l 

1  +  I  1  -  Px  (A)  -  N  +  1  -  l  P  (A) 

j«2  L  j-2  J  j-2  3-2 


N 


N  +  l  -  l  X 

J-l  3 


If  r  simulation  runs  are  made,  we  propose  to  estimate  E[N]  by  the  average 
of  the  above  single  run  estimators. 

As 


E [E [ I j  |  X0,  ....  Xj_2]]  ■=  Ellj]  , 


it  follows  that  our  estimator  is  unbiased.  Also  by  the  conditional  variance 
formula 


Var  (EUj  |  X0,  ....  X^J)  <  Var  (1^) 


However,  the  above  does  not  imply  that 


Var 


(X  EIIJ  1  X0 . Xj-2])  i  *“  (j2  Jj) 


We  now  prove  the  above  inequality,  thus  showing  that  our  estimate  has  smaller 
variance  than  N  . 


Proposition  1; 


Var 


E[Ij  |  XQ,  ...,  Xj  2lj  <  Var  [N]  . 


Proof: 


Let  •  {Xq,  ...»  X±}  denote  the  history  up  to  time  i  and  define 


*  r .  ***  •  '  v; 

O  •’*  •**  /*.A  , 

i*  V-  v'  v  jwW 


h  "  E[Ii  ‘  Hi-2]  * 


m  •  »  lj*  ■»  «  ■»  ■■  m,  a  »  «  ■ ’  m  *  •  •  •  * 

•-  •  k.-. A*-'- 

*•  t  \  fc*'  **'  ***  ^  •*'  "  »  *  ***  ***  .  *.1' 


PMS-JJ 


V  , 

1  '*  • 


•  V’i' 

*>  <*• _ _ 

r~m 


I 

•  ••'J 


[]  y  ■  \ vat  <v +  y  w 


Var  jj  l\  =  l  Var  (ij)  +  i[  Cov  . 

As  Var  (!L)  £  Var  (1^)  and  E [ I ^ ]  =  E [ X ^ ]  ,  it  suffices  to  prove  that 

<  E[I±I  ]  for  all  i  +  j  . 


However,  for  i  <  j  , 


eU^] 


EIEtI1  I  I 

EtEiiilj  I  Hj.j]] 

EliiV 


E'V 

E|IiV 


since  0  <_  1.^  <_  1 


since  1.1,  ■  1^ 
i  J  J 


which  proves  the  result. 


Remark: 


It  is  illuminating  to  write  the  proposed  estimator  of  E(N)  as 


As  £  A  ,  the  total  hazard,  has  mean  1, the  above  is  (as  we  know)  unbiased. 

j=l  j 

Intuitively,  it  has  smaller  variance  than  N  for  in  simulation  runs  for  which 
N  is  small  (large)  the  total  hazard  will  also  (usually)  tend  to  be  small 
(large)  and  so  the  additional  term  in  parenthesis  acts  as  a  correction  factor. 
Of  course,  an  even  better  estimator  is 


" + 4  •  l 


where  a 


Li  j 
■r  0  \) 


,  but  the  difficulty  with  this  is,  of  course. 


that  a  is  usually  not  known. 


1.2  Method  2 


Suppose  that  we  make  a  total  of  r  simulation  runs.  For  run  i  of 


length  ,  define  for  n  <_ 


Xn  *  P  i  (A> 

n  x1  , 

n-1 


where  Xn_^  is  the  (n  -  l)st  state  visited  in  the  ith  run.  Now,  for  a  given 
n  ,  consider  each  run  in  which  >  n  and  let 


i:N(i) 


I  A1/, 

(i)._  n/ 


where  r(n)  «*  number  of  i 


:  N(i)  >  n  . 


Hence,  A  is  an  estimate  of  the 
n 


discrete  hazard  (or  failure)  rate  value  P{N  ■  n  N  >  n)  .  As 


E[N]  =  I  P(N>  j}+l 
j=2 

®  j_x 

=  l  n  P{N  >  n  |  N  >_  n}  +  1  , 
j*2  n=l 


we  have  as  a  point  estimate  of  E[N] 


oo  j-1 

y  n  (1  -  I  )  +  1 

•*  -  n 

j*2  n=l 


where  A  is  taken  to  equal  1  tor  n  >  max  N 
n  i 


1.3  Numerical  Example 


Let  {X^  ,  n  >_  0}  be  a  Markov  chain  having  the  following  transition 


matrix: 


0 

1 

2 

3 

f.  23 

.5 

.25 

.02  > 

.2 

.4867 

.3 

.0133 

.6 

.2 

.19 

.01 

Vo 

0 

0 

1  ) 

Letting  Xq  *  0  ,  the  process  {Xr}  has  been  simulated  1000  times  to  deter 
mine  the  expected  first  passage  time  to  state  4. 


standard 

deviation 

E[N]  of  estimate 


68.281 

66.988 

67.015 

68.059 


2.1  Method  1 


As  shown  in  Section  1.2,  the  quantity  X^  is  an  estimate  of 
,  j 

P{N  ■  n  |  N  ^  n}  ,  and  thus  H  (1  -  X  )  is  an  estimator  of  P{N  >  j}  . 

1X1 

To  approximate  the  variance  of  this  estimator,  we  will  make  use  of  the 

so-called  Delta  Method  used  to  approximate  the  mean  and  variance  of  a  func- 

2 

tion  g(X)  when  E[X]  =  p  ,  Var  (X)  =  a  are  known.  This  method  approx¬ 
imates  g(X)  by 

g(X)  «  g(p)  +  g'(p)(X  -  p) 


leading  to 


E(g(X))  a  g(p) 


Var  (g(X))  a  (g’(u))2a2 


We  now  show  how  to  approximate  Var 
Let 


(i a  - *■>) 


k 

p  =  l  -  x  ,  m.  -  n  p  . 

n  n  k  .  n 

n=l 

That  is,  to  obtain  Pr  ,  one  considers  each  of  the  r(n)  simulation  runs 

which  do  not  end  before  n  transitions  and  for  the  runs  it  takes  the  average 

probability  that  the  nth  transition  does  not  take  the  process  into  A  . 

Hence,  when  defined,  P  is  an  unbiased  estimate  of  P  =  P{N  >  n  I  N  >  n}  . 

n  n  1  — 


log  Mj^  -  I  log  P 
n=l 


and  so 


Var  (log  M.  )  a  £  Var  (log  P  ) 
n-1  n 


a  l  Var  (P„)/P„ 
n-1  n  n 


where  the  first  approximation  results  from  assuming  that  the  P 
independent  and  the  second  from  applying  the  Delta  Method.  As 


Var  (P  )  -  Var  (1  -  P  ) 
n  n 


E[PX  (A>  I  «  >  n  -  lj  -  q*  /r(n) 


where 


qn  -  1  -  PR  -  P{N  ■  n  |  N  >  n}  ,  we  have 


EP:  (A)  I  N  >  n  - 


(log  M^)  a 


n-1 


p2r(n) 


Again,  using  the  Delta  Method, we  have 


Var  (log  Mk)  a  Var  (M^/E  tMR] 


and  so 


.*•  /•  .  *•  ;*•*  /•  .*•  »v  .*■  /*  >  -1'-  .*•  .  .  •  .'* 


j 


(A)  N  >  n 


Hence,  we  can  obtain  a  rough  estimate  of  the  variance  of  our  estimate 
by  using  the  above  and  approximating  E[M. ]  ,  E  P?  (A)  |  N  >  n  -  lj  , 

j  Ln-l  J 


P  ,  q  by 
n  nn  J 


E[Mj ]  C-t  Mj 


[p*  (A)  |  N  >  n  -  1 

1  6gt  I 

(Pi  (AA 

L  n-1 

J  i:N^>n 

\Xn-l  / 

l  P,  (A)  /  r(n) 
i:N  >n  X  .  ' 


i-  n-1 


q  C-t  1  -  P 
Mn  r 


where  ■  means  "is  estimated  by". 

2.2  Method  2 

A  second  method  for  estimating  P(N  j)  is 


t  l  Aj(i) 

i-1  J 


where  A 


«  -  {  , 
J  nil  n 


and  where  X^  ■  0  if  <  n  . 

n 


It  is  unbiased  and  since  the  different  simulation  runs  are  independent 


its  variance  is  —  Var 


U.»\  . 


Proposition; 


Var 


m  2E|Aj(1)l{N(1)  <  j}  -EX  (1)1{N(1)  <  j} 


-  {P(N(1)  <  j)}2 


where  !{•}  is  the  indicator  function  of  the  specified  event  {•} 


Proof : 


(JL  ‘i")' 

•  •  i  k  «”  ■  k  Hi 

-2  l  A(1)X(1)  -  l  (x(1))  . 

i  n  n  ‘‘i  \  n  / 

n-1  n-1  '  ' 


Taking  expectation  of  both  sides,  we  obtain 


[(A™)2]  * 2 1  EK1>i“>]  - 1  E[(X”U)2] 


Using  the  basic  property  of  conditional  expectation,  we  get 


[(Aji>)2]  ■ 2  X  eKi>i{n<i>  •  ">]  -  j,  e[>»i>iik 

■  2E[A‘u>1{K<1>  -  E[^(i)ltK<1>  - Jl] 


(1) 


and  the  result  follows. 


The  above  proposition  provides  good  insight  to  the  condition  under  which 


the  proposed  method  is  superior  to  the  classical  crude  method.  Consider  a 
chain  in  which  the  set  A  is  "rare",  that  is,  for  x  t  A  ,  p  <_  P  (A)  <_  p 
where  0  <  p^  <  are  "small",  then  for  a  moderate  j  ,  P(N^  <_  j)  =  p  is 
small  and  the  variance  which  appears  in  the  crude  method  is  (roughly)  pro¬ 
portional  to  p/r  whereas  the  variance  of  the  proposed  estimator  is  of  the 

2 

order  of  magnitude  of  p  /r  .  For  small  p  ,  this  is  a  significant  improve¬ 
ment. 


2.3  Numerical  Example 

{X  ,  n  >  0}  is  as  in  1.3.  It  was  simulated  1000  to  determine  P(enters 
n  — 

state  3  by  time  10). 


Probability  and  Standard  Deviation 
of  Estimate 


Probability 

Standard 

Deviation 

theoretical 

.1408 

— 

crude 

.1480 

1123  x  10“5 

method  1 

.1407 

5.5  x  10-5 

method  2 

.1408 

93.5  x  io-5 

Relative  efficiencies: 


method  1  to  crude:  205 
method  2  to  crude:  12 
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3.  ESTIMATION  OF  P{XJJ  e  B} 

Let  B  C  A  denote  a  subset  of  A  and  consider  the  problem  of  estimating 

the  probability  that  the  final  state  is  in  B  .  To  do  this,  let  N  equal  the 

number  of  transitions  until  B  is  reached  (and  so  Ng  -  «  if  i  B) . 

Our  approach  to  estimating  P{ e  B}  will  consist  of  working  with  the 

following  modified  version  of  the  hazard  rate  function  of  N — namely, 

B 

A*(n)  -  P{Nb  -  n  |  N  >  n}  . 


To  estimate  A*(n)  ,  consider  those  simulation  runs  for  which  >  n  -  1 

B 

and  then  take  as  the  estimate  of  A*(n)  the  average  value  of  Pv  (B)  for 

n-1 

those  runs.  That  is,  AD(n)  ,  the  estimate  of  A*(n)  is  given  by 

B  B 


A*(n) 


:N(i*>n 


'r  (n) 


where  r(n)  is  the  number  of  runs  for  which  N  >  n  .  Take  A*(n)  to  equal 
0  if  r  (n)  -  0  . 

Hence,  we  have  an  estimate  for  P{Ng  «*  n}/P{N  >_  n}  .  Also,  by  using  the 

method  of  Section  2.1, we  can  estimate  P{N  >_  n}  and  so  by  using  the  ratio 

of  these  estimators  we  have  a  method  for  estimating  P{Ng  -  n}  .  By  summing 

00 

over  all  n  ,  we  can  then  estimate  \  P{Ng  ■  n}  -  P{X^  c  B}  .  In  addition, 

n-1 

as  we  are  assuming,  but  not  using  in  our  estimation  procedure,  the  result 
that  P{X^  e  A}  -  1  ,  we  can  improve  our  estimate  of  P{X^  e  B) — call  it 
P{Xjj  e  B} — by  normalizing.  That  is,  our  estimate  of  P{X^  e  B}  is 
§{3^  c  B}/P{Xjj  e  A}  . 
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Remark: 

The  reason  we  use  the  modified  hazard  rate  function  P{Ng  =  n  |  N  ^  n} 
instead  of  the  usual  hazard  rate  function  P{N„  =  n  I  N.  >  n)  is  that  a 

D  D  — 

direct  estimate  of  this  latter  expression  does  not  use  all  available  informa¬ 
tion.  For  if  Ng  ^  n  then  Px  (B)  estimates  P{Ng  -  n  |  Ng  >_  n}  ,  but  in 

n— 1 

addition  P^  (A  -  B)  yields  valuable  information  about  P{X^  i  B}  .  This 
n-1 

information  is  captured  in  our  method  by  coming  into  play  in  our  estimate  of 
P{N  n}  .  Hence,  our  ratio  estimate  makes  use  at  each  stage  of  both  the 
probability  that  the  next  state  will  be  in  B  and  the  probability  that  B 
will  never  be  reached.  A  product  form  estimate  of  P{N„  <  <*>}  which  uses 

D 

only  direct  estimates  of  the  (unmodified)  hazard  rate  function  will  not 
utilize  this  latter  information. 


.  •.  V  •• 
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4.  CONTINUOUS  TIME  VERSION 

Suppose  now  that  the  Markov  process  is  in  continuous  time-spending  and 

exponential  amount  of  time  with  mean  y(x)  ,  0  <  y(x)  <  ®  ,  when  in  state 

.  -* 

x  before  transiting  into  another  state  according  to  P^  .  Suppose  also 
that  the  Markov  chain  is  regular  in  the  sense  that  with  Probability  1  it 
only  makes  a  finite  number  of  transitions  in  any  finite  time  interval.  Let 
T  denote  the  time  it  takes  for  the  process  to  enter  A  . 

4.1  Estimating  E[T] 

In  analogy  with  Proposition  1, for  any  given  simulation  run,  an  estimate 
with  smaller  variance  and  the  same  mean  as  T  is 

T 

s^(A)/y(X(s))ds 

0 

where  X(s)  is  the  state  at  time  s  .  Now,  let  X^.X^,  . . . ,  Xjj  denote  the 
sequence  of  states  visited.  The  above  estimator  can  be  improved  thusly: 

E[T  +  1  -  total  hazard  |  N.Xq.X.^  ...,  Xjj] 

-  E[T  |  N,Xq,  ....  Xj,]  +  1  -  E[hazard  |  N,XQ,  ....  Xjj] 

[N-l 

-  p (XQ)  +  ...  +  +  1  -  El  l  TtP  (A)/y(X 

N-l 

-  y(Xn)  +  ...  +  yCX^  .)  +  1  -  l  PY  (A) 

i-0  Ai 

where  in  the  above  represented  the  time  spent  in  state  XA  before  the 

transition  to  X^+^  occurs  (and  is  thus  exponential  with  mean  y(X^)  in¬ 


T  +  1  -  total  hazard  -  T  +  1 


- 


X( 


dependent  of  the  next  state  visited) . 


Hence, our  estimate  for  E[T]  in  any  run  is  just  the  conditional  value 
of  T  given  the  sequence  of  states  visited  plus  1  minus  the  total  discrete 


hazard  given  these  states.  The  overall  estimate  is  then  obtained  by  averaging 
this  quantity  over  all  simulation  runs. 

A. 2  Estimating  P{T  <.'  t) 

Again  we  recommend  conditioning  on  the  sequence  of  states  visited 
H,Xq,X^,  ...»  X^  .  That  is,  for  a  given  simulation  run,  let 

r.j1  lf  Tlt 

(0  otherwise  . 

Then, 

E[I  |  N,Xq,  ...,  XN]  -  P{tq  +  ...  +  tn_1  <  t} 

where  Tq,  ...»  TN_^  are  independent  exponentials  with  respective  means 
m(Xq)  ,y(X^) ,  ...,  p(Xjj_^)  .  Where  appropriate  (N  relatively  large  and  none 
of  the  m(Xq),  ...,  y(X^_^)  dominant)  the  above  can  be  approximated  by  the 
central  limit  theorem  and  where  this  is  not  appropriate  a  simulation  can  be 
done  using  the  approach  of  Section  A. 2.1. 

A. 2.1  Simulation  Estimation  of  a  Convolution 

We  start  with  a  proposition  which  will  be  necessary  in  the  sequel. 

Proposition: 

Let  X  ■  (X^,X2»  ...)  be  any  sequence  of  random  variables  and  let 

gi(X)  denote  arbitrary  functions,  i  ■  1,  ...»  n  .  Then,  for  >_  0  , 
n 

7  X .  ■  1  ,  we  have  that  for  any  random  variable  Z 


Var  \k  X±ElZ  ^  8i-)]j  1  VSr  (Z)  * 


Proof; 

Let  M  denote  a  random  variable  independent  of  Z  and  of  X  and 
such  that 


P{M  *  i}  =  X±  ,  i  =  1 . n  . 


Then,  by  the  conditional  variance  formula, 

Var  (Z)  >  Var  (E[Z  |  M.g^X)]) 

-  E[Var  (E[Z  |  M.g^X)]  |  X)] 

+  Var  (E[E[Z  |  M.g^X)]  |  X] ) 
>  Var  (E[E[Z  |  M.g^X)]  |  X]) 

-  Var  V12  I  8i(Pl)  -M 


Now,  let  X^,X2»  . . . ,  Xn  be  independent  with  X^  having  distribution 

F.  and  suppose  we  want  to  use  simulation  to  estimate  P{X,  +  ...  +  X  <  t} 
i  In  — 

To  do  so,  let 


Z  -  I 


1  if  l  X  <  t 

i-1  1 

0  otherwise  . 


We  can  use  the  above  by  letting 


g ±<x)  -  l  X 
i+i  J 


and  so 
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EH  I  8l(X)l  -  Ft(t  *  I  Xj) 


and  so  an  unbiased  estimate  of  the  convolution  having  smaller  variance  than 


I  is  given  by 


t,  +  ...  +  Xn  <_  t>  e«t  l  >1Fi(t  •  L  Xj) 


where  A  >_  0  ,  £  A  =  1  . 

1  1=1  1 

As  the  above  estimator  is  an  improvement  over  I  for  all  nonnegative 

vectors  _A  that  sum  to  unity,  it  remains  to  present  a  method  for  picking 

the  "optimal"  A  .  To  do  so,  let  us  recall  that  if  W,,W_,  ...»  W  are 
—  l  l  n 


independent  random  variables  having  common  mean  then  variance  of  £  A.^ 

1/Var  (Wi)  /  1“1  v 

is  minimized  by  letting  A  ■  -  .  Whereas  F.  (t  -  \  X  I  , 

l  l/Var  (W.)  *'  J"1  S' 

j=l  J 

1  =  1,  ...,  n  are  clearly  not  independent,  let  us  nevertheless  use  this 
result  to  choose  the  convex  linear  coefficients  A^  .  Hence,  we  must 


determine  Var  F 


‘(c  - 1 x;) 


Now,  if  Y  is  a  random  variable  and  if  F  is  a  distribution  function 
with  density  F'  =  f  ,  then  by  the  Delta  Method 


Var  (F(t  -  Y))  2S  Var  (Y)[f(t  -  E[Y])]  . 


Hence , 


Var 


M'  '  jh  Xj)]  =  [lk  '  ill  Uj)]2  sli  t 


We  thus 


where  y^  and 
recommend  taking 


are  respectively  the  mean  and  variance  of 


l  '*) 1 

i  n 

/  2/ 

v  ?\-l 

I 

fi(t  -  l  v. 

1  s  °i 

i*l 

\  \  1H  ' 

7  j*i  j/ 

In  the  special  case  where  is  exponential  with  mean  y^  and  so 


V* 


1  “x/uj  2  2 

—  e  ,  o.  -  p,  ,  we  have 

Mj  J  i 


where  C  is  chosen  to  make 


A. 2. 2  Estimating  P{T  <.  t} 


n 


by  Uniformizing 


Another  possibility  for  estimating  P{T  <_  t}  is  to  uniformi2e  the 

Markov  process — which  can  be  done  assuming  that  inf  y(x)  >  0  .  This  trans- 

x 

forms  the  Markov  process  into  one  which  spends  an  exponential  time  with  rate 
y(x)  =  y  =  1/X  in  state  x  .  (See,  for  instance,  Ross  [3]).  Hence,  we  can 
compute  P{T  £  t}  by  conditioning  on  the  number  of  transitions  it  takes  to 
enter  A  .  This  gives 
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P{T  <_  t}  =  l  P{T  <_  t  |  N  -  n}P{N  *  n} 
n 


-  I  G  (t)P{N  -  n} 
«  n,A 


where 


i 


Ae 


-Ax  (Ax) 


n-1 


(n  -  1)  ! 


dx 


1  - 


-At 

e 


n-1 


I 

j-o 


(At)j 

j! 


is  the  gamma  (n,A)  distribution  function.  Hence,  we  can  estimate  P{T  <_  t) 
by  uniformizing  the  process  and  then  use  the  resultant  simulation  runs  to 
estimate  P{N  =  n}  .  This  last  quantity  can  be  estimated  either  by 


P{N 


if  N  >_  n 
if  N  <  n 


and  then  averaging  over  all  runs,  or  by  using 


P{N 


n-1 

n)  e®  A  n 
i-1 


n 


(1  -  A±)  . 


Remark: 

By  "uniformizing",  we  eliminate  all  variance  due  to  the  times  spent  in 
each  state.  However,  we  pay  the  price  of  increasing  the  variance  on  the 
number  of  transitions  involved.  This  procedure  is  thus  only  recommended 
when  the  values  p(x)  are  relatively  constant — for  in  this  case  very  few 
additional  transitions  result. 
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